Simulation and detection of Dirac fermions with cold atoms in an optical lattice 
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We propose an experimental scheme to simulate and observe relativistic Dirac fermions with cold 
atoms in a hexagonal optical lattice. By controlling the lattice anisotropy, one can realize both 
massive and massless Dirac fermions and observe the phase transition between them. Through 
explicit calculations, we show that both the Bragg spectroscopy and the atomic density profile in a 
trap can be used to demonstrate the Dirac fermions and the associated phase transition. 
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Control of ultracold atoms in an optical lattice opens 
up many avenues to explore some fundamental phenom- 
ena at the forefront of condensed matter physics [l], [I, S 
0, 5]. By designing configurations of this atomic system, 
one can simulate effective theories that are very differ- 
ent from the microscopic atomic physics. In this paper, 
we add an unusual example to the avenues of quantum 
simulation by showing that ultracold atoms in an optical 
lattice can be used to investigate physics associated with 
relativistic Dirac fermions. The ultracold atomic gas, 
as the coldest setup in the universe, is one of the most 
non-relativistic systems. Nevertheless, we will see that 
effective theories for the quasiparticles in this system can 
become relativistic under certain conditions. 

We simulate Dirac fermions with single-component 
cold atoms in a two-dimensional hexagonal lattice. This 
lattice can be formed through interference of three laser 
beams, as we show below. The physics here is closely 
related to the properties of electrons in the graphene ma- 
terial formed with a single layer of carbon atoms 0, 



B & 11, m. The graphene, with its emergent rela- 
tivistic massless quasiparticles, has recently raised strong 
interest in condensed- matter physics 0, Isl, [ul, [l2| • 
Compared with the graphene, the system with cold atoms 
in an optical lattice may offer more controllability. For 
instance, we show that one can realize both massive 
and massless Dirac fermions by controlling anisotropy of 
the optical lattice. This anisotropy can be conveniently 
tuned through variation of the trapping laser intensity. 
Under such a tuning, one can also observe a quantum 
phase transition in this system. This phase transition is 
not associated with any usual symmetry breaking, but 
instead it is characterized by a topological change of the 
fermi surfaces [HI, [l3]- To detect the massive and the 
massless Dirac fermions and the phase transition between 
them, we calculate the Bragg spectrum for this system 
as well as its atomic density profile in a trap. From this 
calculation, we show that the conventional atomic detec- 
tion techniques based on the Bragg spectroscopy [l5[ or 
the density profile measurement [1, 16', jl^ can be used to 
demonstrate the Dirac fermions and the phase transition. 
For cold atoms, one realizes an effectively two- 



dimensional system by raising the potential barrier of 
the optical lattice along the z direction to suppress the 
vertical tunneling between different planes. Then, in the 
x-y plane, one can form a hexagonal optical lattice with 
three laser beams. For instance, as shown in Ref. 0], 
with three detuned standing- wave laser, the optical po- 
tential is given by 

cos6>j+?/sin6>j)+7r/2] (1) 

where Oi = 7r/3, 62 = 27r/3, 0^ = 0, and kL is the optical 
wave vector. If Vi = V2 = V^^ the energy contour of the 
potential V{x^ y) is shown in Fig. 1(a), where its minima 
(marked with the solid dots) form a standard hexago- 
nal lattice. It is easy to tune the potential barriers Vj by 
variation of the laser intensities along different directions. 
With different Vj , one can still get a hexagonal lattice but 
with a finite anisotropy. For instance. Figure 1(b) shows 
the potential contour with Vi = V2 = O.9IF3, where the 
neighboring sites along the horizontal direction have a 
larger distance and a higher barrier. As the atomic tun- 
neling rate in an optical lattice is exponentially sensitive 
to the potential barrier, this control provides an effective 
method to tune the anisotropy of the atomic tunneling 
rate in this lattice. 

A hexagonal lattice consists of two sub-lattices denoted 
by A and B as shown in Fig. 1(c). We consider single- 
component fermionic atoms (e.g., ^^K^^Li etc.) in this 
hexagonal lattice. For single-component fermions, the 
atomic collisions are negligible at low temperature. The 
Hamiltonian is given by the simple form 



H = 



{albj^h.c), 



(2) 



where (i, j) represents the neighboring sites, and and 
bj denote the fermionic mode operators for the sub- 
lattices A and B, respectively. The tunneling rates 
tij in general depend on the tunneling directions in an 
anisotropic hexagonal lattice, and we denote them as 
^i,^2,fe corresponding to the three different directions 
(see Fig. 1(c)). In this work, for simplicity, we assume 
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ti = t2 = t and ts = /3t, where (3 > is the anisotropy 
parameter. In the fohowing, we show that by tuning the 
anisotropy /3, the quasiparticles in this system change 
their behaviors from massless to massive Dirac fermions, 
with a quantum phase transition between the two cases. 

The Hamiltonian (2) can be diagonahzed with a simple 
extension of the method for the graphene material [l2|. 
For the sub-lattice A, the positions of the sites can be 
expressed as A = miai + m2a2 , where mi and m2 are 
integers, and the basis vectors ai = (V^, ~^){^/'^)^ ^2 = 
(0,a) (a = 27r/ {V^kL) is the lattice spacing). The sites 
in the sub-lattice B are generated by a shift B = A + b 
with three possible shift vectors bi = (l/>/3, l)(a/2), 
b2 = (l/V3,-l)(a/2), and bs = {-a/V3,0) (see Fig. 
1(c)). The first Brillouin zone of this system also has a 
hexagonal shape in the momentum space with opposite 
sides identified but rotated an angle of 7r/6 relative to 
the hexagon of the real-space lattice (see Fig. 1(d)). Cor- 
responding to two different sites A and B in each cell 
in real hexagonal lattice, only two of the six corners in 
Fig. 1(d) are inequivalent, usually denoted as K and K' . 
One can choose K = (27r/a)(l/v^, 1) and = -K. 
With a Fourier transform a] = --^^j^exp(zk • A^)aj^ 

and b^j = ^j^exp(zk • Bj)6j^, where N is the num- 
ber of sites of the sub-lattice A (or 5), the Hamilto- 
nian (2) simplifies to = X]k[^(^)^k^k + h.c], where 
0(k) = -EU^tsexp{ik-hs). The energy eigenvalues of 
H are given by = ±|^(k)|, which has the expression 
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FIG. 1: (Color online) The hexagonal lattices, (a) and (b) 
The contours with three potentials described in Eq.([T]). The 
minima of the potentials are denoted by the solid dots. All 
are the same in (a) and = V2 = 0.91^3° in (b). (c) Decom- 
position of the hexagonal lattice as two triangular sublattices 
A and B with anisotropic tunnelling, (d) The Brillouin zone 
of the hexagonal lattice. The dispersion relations are shown 
in (e) for (3—1 (gap less state) and (f) for [3 — 2.^ (gapped 
state) . 



Qx and Qy^ the dispersion relation ([3|) becomes 



2 cos{kya) + 4(3 cos (^V3kxa/2^ cos {k^ 

(3) 



The energy versus the momentum from Eq.Q is plot- 
ted in Fig. 1(e) and 1(f) for different values of the 
anisotropy (3. There are two branches of curves corre- 
sponding to the ± sign in Eq. (3). When < < 2, the 
two branches touch each other, and around the touching 
points there appears a Dirac cone structure. One has the 
standard cone same as the graphene material with (3 = 1 
S S [lH 5 and as p deviates from 1 , the cones squeeze in 
X or y direction, but they still touch each other. When 

> 2, the two branches separate with a finite energy 
gap A = \t\{P — 2) between them. So, across the point 
(3 = 2^ the topology of the fermi surface changes and 
there is a corresponding quantum phase transition, al- 
beit no symmetry breaks at this point. With this phase 
transition, the system changes its behavior from a semi- 
metal to an insulator at the half filling case (half fill- 
ing means one atom per cell; note that each cell has 
two sites). Around the half filling, the fermi surface is 
close to the touching points and one can expand the mo- 
mentum k around one of the touching points {k^^ky) as 
{kx, ky) = {k^ ^ qx^ky ^ qy). Up to the second order of 



(4) 
a/2). 

where A = 0, ^/^(3ta/2 and ta^l - 

for < /3 < 2; and A = \t\ {(3 - 2), = tay/3(3/2 
and Vy = tay^ (3/2 — 1 lor (3 > 2. This simplified dis- 
persion relation is actually a good approximation as 
long as qx^qy ^ l/2a. We see that represents the 
standard energy-momentum relation for the relativistic 
Dirac particles, with A taking the meaning of mass, and 
and Vy replacing the light velocity. The wavefunction 
for the quasiparticles around the half filling then satisfies 
the Dirac equation ihdt^ = H^^^ where the relativistic 
Hamiltonian Hd is given by 



Hd = aoA + v^axPx + v^ayPy, 



(5) 



where the (/i = 0, x, y) matrixes satisfy the Grassma- 
nian algebra a^a^f + (^fi'(^fi = 2(5^^/, and for the 2+1 
dimensional system, they can be taken as the three Pauli 
matrixes (Jz^cTx^cry [12]. 

In the above, through an analogy to the graphene 
physics, we have shown how to realize massive and mass- 
less Dirac fermions with cold atoms in an anisotropic 
hexagonal lattice. A more important question for this 
system, however, is how to experimentally verify the ex- 
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istence of these relativistic quasiparticles and the asso- 
ciated quantum phase transition. The detection method 
for ultracold atoms is very different from that for con- 
densed matter materials, and the widely used technique 
based on the transport measurements for the latter is typ- 
ically not available for the atoms. Nevertheless, there are 
some specific detection methods for the trapped atomic 
gas, and in the following we show how to confirm the 
relativistic quasiparticles and the phase transition with 
the density profile measurement [l|, ll6|, ll7[ and the Bragg 
spectroscopy [15j. 

The density profile of the trapped atoms can be mea- 
sured through the time-of-flight imaging with the light 
absorption. For free fermions, we have ballistic expan- 
sion, and from the final measured absorption images, one 
can reconstruct the initial real-space density profile of the 
trapped gas [17]. Now we show this density profile pro- 
vides critical information for both massive and massless 
Dirac fermions. For trapped fermions, the local density 
approximation (LDA) is typically well satisfied, and un- 
der the LDA, the local chemical potential varies with the 
radial coordinate by /i = /iQ — V'(r), where /io is the chem- 
ical potential at the trap center and V(t) = mcc;^r^/2 is 
the global harmonic trapping potential [l8|. So, /i is a 
monotonic function of r, and the density profile n(r) is 
uniquely determined by the equation of the state n{fi). 

For this system at temperature T, the atomic density 
(the number per unit cell) is given by 



li^i) = j-l f{k. 



(6) 



where So = Sir'^/VSa'^ is the area of the first Bril- 
louin zone of the honeycomb lattice, and f{kx^ky^ ji) = 
l/{exp[(£^k — + is the Fermi distribution. At low 
temperature (T ~ 0), this density profile n(/i) is shown 
in Fig. 2(a) and 2(b) for the parameters with massless 
and massive Dirac quasiparticles, respectively. One can 
clearly see that for the gapped phase with massive Dirac 
fermions, one has a plateau at the atom density n = 1 
in the density profile. For the case with massless Dirac 
fermions, there is no such a plateau. So the plateau is as- 
sociated with massive quasiparticles, and its emergence 
with tuning of the lattice anisotropy provides an unam- 
biguous signal for the quantum phase transition between 
the two cases. 

To further confirm the massless Dirac fermions, one 
needs to have evidence for their linear dispersion rela- 
tion with the Dirac cone structure. At T ~ 0, n(/i) 
in Eq. (6) simplifies to n(/i) = (l/5'o) <^^dkxdky. 
Around the half filling (n = 1 which corresponds to 
the touching point in Fig. 1(e)), with a variation Sjii 
in the chemical potential, ^n(/i) geometrically represents 
the cross section of the Dirac cone, so it must be pro- 
portional to ((5/i)^. In Fig. 2(c) and 2(d), we show the 
numerically calculated derivative of n(/i), and indeed at 



the vicinity of the half filling, we see that ^ is linearly 
proportional to S/j. with the explicit asymptotic expres- 
sion ^ = where we have assumed the velocity 
'^x = = ^- So, experimentally, from the measured den- 
sity profile n(/i), one can determine its slope. The latter 
quantity, with its linear form shown in Fig. 2(d), sig- 
nals the linear dispersion relation around the Dirac cone, 
which confirms the massless Dirac fermions. 
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FIG. 2: (Color online) The number density of atoms n per 
unit cell as a function of the chemical potential /i (correspond- 
ing to a rescaled atomic density profile in a trap) for (a) = 1, 
and (b) /3 = 2.5. A plateau with a width 2/3 — 4 appears for 
the latter case which corresponds to the case when the chemi- 
cal potential sweeps inside the energy gap. (c) The derivative 
dn/d/j. as a function of the chemical potential /i for /3 = 1. 
(d) An enlarged part of dn/d/j at the vicinity of /i = 0. The 
linearity of the curve shows the linear dispersion relation for 
the quasiparticles. 

The Bragg spectroscopy can provide an alternative and 
complementary method to confirm the linear dispersion 
relation for the massless Dirac fermions and the energy 
gap for the massive ones. In Bragg spectroscopy [IB], one 
shines two laser beams on the atomic gas as shown in 
Fig. 3a. By fixing the angle between the two beams 
(thus fixing the relative momentum transfer q = k2 — ki , 
where denotes the wave vector of each laser beam), 
one can measure the atomic (or photonic) transition rate 
by scanning the laser frequency difference uj = UJ2 — ^i- 
From the fermi's golden rule, this transition rate basically 
measures the following dynamical structure factor [lH 

5(q,c^)= ^ |(/k,|i/BKk,)|'^[^c^-^/k.+^.kJ, (7) 

ki,k2 

where Hb = Eki,k2 ^e^^ ''Kki)(/k2 1 + h.c. is the light- 
atom interaction Hamiltonian, and |zki) and l/ks) denote 
the initial and the final atomic states with the energies 
£^iki and Ef\^^ and the momenta ki and k2, respectively. 
At the half filling, the valence band (the lower half of Fig. 
(le) and 1(f)) is fully occupied, and the conduction band 
(the upper half) is empty. In that case, the excitations 
are dominantly around the touching point, and we can 
use the approximate dispersion relation in Eq.(|3j). For 
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the isotropic case (/? = 1) with massless Dirac Fermions, 
we find S{q^u) has the expression 



r 0, 



{(jO < UOr) 
{u > UJr) 



(8) 



where cOr = qv/h {q = |q|) and qr = hco/v. This dy- 
namical structure factor is shown in Fig. 3(b). Note 
that in this case, the lower cutoff frequency cOr is lin- 
early proportional to the momentum difference q, and 
LJr vanishes when q tends to zero. The ratio between 
UJr and q gives the fermi velocity an important pa- 
rameter as the analogy of the light velocity for conven- 
tional relativistic particles. For the anisotropic case with 
> 2, the spectrum in Eq. (|4]) becomes quadratic with 
E « ±{A-^h'^ql/2mx-^h'^qy/2my) for smah momentum 
transfer q, where the effective mass rrix^y = h'^A/v'^^^. 
The dynamic structure factor in this non-relativistic limit 
becomes 



quantum phase transition with cold atoms in a (generally 
anisotropic) hexagonal optical lattice. The characteris- 
tic dispersion relations for the massless or the massive 
Dirac fermions can be confirmed through either the den- 
sity profile measurement or the Bragg spectroscopy. The 
phase transition can be identified with the appearance 
of a plateau in the density profile as one tunes the lat- 
tice anisotropy. The appearance of relativistic quasipar- 
ticles, together with control of gauge fields for this system 



19l . |20| , opens up the prospect to use the ultracold atoms 



to simulate some high energy physics. 



This work was supported by the NSF under grant num- 
ber 0431476, the ARDA under ARO contracts, the A. P. 
Sloan Foundation, the State Key Program for Basic Re- 
search of China (No. 2006CB921800), NCET and NSFC 
(Nos. 10674049, 10474034, and No. 60421003). 



f 0, {u< u^^y) 

S{q.u) = ^ ^ (c^>c.-'^) (9) 

where cj^'^ = 2 A + h?q^ y/Amx^y Its form is shown in 
Fig. 3(b). The lower cutoff frequency uo^'^ in this case 
does not vanish as the momentum transfer go to zero. 
This distinctive difference between the dynamical struc- 
ture factors in Eqs. ([8]) and ([9]) can be used to distinguish 
the cases with massive or massless Dirac fermions. From 
the variation of the cutoff frequency cj^'^ as a function of 
the momentum transfer qx,y^ one can also experimentally 
figure out the important parameters such as the energy 
gap A and the effective masses rrix and rriy. 




FIG. 3: (Color online) Schematic of the Bragg scattering, 
(a) Atoms are illuminated by two laser beams with wave vec- 
tors ki, k2 and frequencies cji, a;2, respectively, (b) The dy- 
namic structure factors S{q^ uj) (not scaled) for the massless 
Dirac fermions (solid line) and for the massive ones in the 
non-relativistic limit (dotted line). The experimentally mea- 
surable quantities uOr and a;^'^ give important parameters for 
the quasiparticles. 

In summary, we have proposed an experimental scheme 
to simulate and observe relativistic Dirac fermions and a 
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